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This is a review of recent advances in our understanding of how Andreev reflection at a super- 
conductor modifies the excitation spectrum of a quantum dot. The emphasis is on two-dimensional 
impurity-free structures in which the classical dynamics is chaotic. Such Andreev billiards differ 
in a fundamental way from their non-superconducting counterparts. Most notably, the difference 
between chaotic and integrable classical dynamics shows up already in the level density, instead 
of only in the level-level correlations. A chaotic billiard has a gap in the spectrum around the 
Fermi energy, while integrable billiards have a linearly vanishing density of states. The excitation 
gap -Bgap corresponds to a time scale h/Ega_p which is classical (?i-independent, equal to the mean 
time Tdweii between Andreev reflections) if Tdweii is sufficiently large. There is a competing quan- 
tum time scale, the Ehrenfest time te, which depends logarithmically on h. Two phenomenological 
theories provide a consistent description of the rg-dependence of the gap, given qualitatively by 
_Egap — min(?i/rdwcii, ^/t"e). The analytical predictions have been tested by computer simulations 
but not yet experimentally. 

PACS numbers: 74.45.-l-c,05.45.Mt,73.23.-b,74.78.Na 



I. INTRODUCTION 

Forty years ago, Andreev discovered a peculiar prop- 
erty of superconducting mirrors _!]■ As illustrated in Fig. 
n an electron that tries to enter a superconductor com- 
ing from the Fermi level of a normal metal is forced to 
retrace its path, as if time is reversed. Also the charge 
of the particle is reversed, since the negatively charged 
electron is converted into a positively charged hole. The 
velocity of a hole is opposite to its momentum, so the 
superconducting mirror conserves the momentum of the 
reflected particle. In contrast, reflection at an ordinary 
mirror (an insulator) conserves charge but not momen- 
tum. The unusual scattering process at the interface be- 
tween a normal metal (N) and a superconductor (S) is 
called Andreev reflection. 

Andreev reflection is the key concept needed to under- 
stand the properties of nanostructures with NS interfaces 
yl . Most of the research has concentrated on transport 
properties of open structures, see Refs. |3il4 ^^r reviews. 
There experiment and theory have reached a compara- 
ble level of maturity. In the present review we focus on 
spectral properties of closed structures, such as the quan- 
tum dot with superconducting contacts shown in Fig. |2] 
The theoretical understanding of these systems, gained 
from the combination of analytical theory and computer 
simulations, has reached the stage that a comprehensive 
review is called for — even though an experimental test 
of the theoretical predictions is still lacking. 

An impurity-free quantum dot in contact with a su- 
perconductor has been called an "Andreev billiard" [3-^ 





^ Open structures containin; 
called "Andreev billiards' 
selves to closed systems. 



an antidot lattice have also been 
, but in this review we restrict our- 



FIG. I: Normal reflection by an insulator (I) versus Andreev 
reflection by a superconductor (S) of an electron excitation in 
a normal metal (N) near the Fermi energy Ef- Normal reflec- 
tion (left) conserves charge but does not conserve momentum. 
Andreev reflection (right) conserves momentum but does not 
conserve charge: The electron (e) is reflected as a hole (h) 
with the same momentum and opposite velocity. The missing 
charge of 2e is absorbed as a Cooper pair by the supercon- 
ducting condensate. The electron-hole symmetry is exact at 
the Fermi level. If the electron is at a finite energy E above 
Ep, then the hole is at an energy E below Ep- The energy 
difference of 2E breaks the electron-hole symmetry. From 
Ref. i. 



The name is appropriate, and we will use it too, because 
it makes a connection with the literature on quantum 
chaos 0, H • A billiard (in the sense of a bounded two- 
dimensional region in which all scattering occurs at the 
boundaries) is the simplest system in which to search for 
quantum mechanical signatures of chaotic classical dy- 
namics. That is the basic theme of the field of quantum 
chaos. By introducing a superconducting segment in the 
boundary of a billiard one has the possibility of unravel- 
ing the chaotic dynamics, so to say by making time flow 
backwards. Andreev billiards therefore reveal features of 
the chaotic dynamics that are obscured in their normal 
(non-superconducting) counterparts. 




FIG. 2; Quantum dot (central square of dimensions 500 nm x 
500 nm) fabricated in a high- mobility InAs/AlSb heterostruc- 
ture and contacted by four superconducting Nb electrodes. 
Device made by A. T. Filip, Groningen University (unpub- 
lished figure). 



The presence of even the smallest superconducting seg- 
ment suppresses the quantum mechanical level density at 
sufficiently low excitation energies. This suppression may 
take the form of an excitation gap, at an energy -Bgap well 
below the gap A in the bulk superconductor (hence the 
name "minigap"). It may also take the form of a level 
density that vanishes smoothly (typically linearly) upon 
approaching the Fermi level, without an actual gap. The 
presence or absence of a gap is a quantum signature of 
chaos. That is a fundamental difference between normal 
billiards and Andreev billiards, since in a normal billiard 
the level density can not distinguish chaotic from inte- 
grable classical dynamics. (It depends only on the area 
of the billiard, not on its shape.) 

A powerful technique to determine the spectrum of a 
chaotic system is random- matrix theory (RMT) |3.l9l[lO|. 
Although the appearance of an excitation gap is a quan- 
tum mechanical effect, the corresponding time scale 
^/Egap as it follows from RMT is a classical (meaning 
?i-independent) quantity: It is the mean time Tdwcii that 
an electron or hole excitation dwells in the billiard be- 
tween two subsequent Andreev reflections. A major de- 
velopment of the last few years has been the discovery of 
a competing quantum mechanical time scale te oc | ln?i|. 
(The subscript E stands for Ehrenfest.) RMT breaks 
down if Te ^ Tdwoii and a new theory is needed to de- 
termine the excitation gap in this regime. Two different 
phenomenological approaches have now reached a con- 
sistent description of the r^j-dependence of the gap, al- 
though some disagreement remains. 

The plan of this review is as follows. The next four sec- 
tions contain background material on Andreev reflection 
(Seem, on t^^s minigap in NS junctions (Sec. IIII|I . on 
the scattering theory of Andreev billiards (Sec. lIV|) . and 
on a stroboscopic model used in computer simulations 
(Sec. E|). The regime of RMT (when te < Tdwoii) is de- 
scribed in Sec. IVII and the quasiclassical regime (when 



Te ^ T-dwcii) is described in Sec. lVIII The crossover from 
Egap — ?i/T'dwcii to Egap ~ ^/te is the topic of Sec. IVIIII 
We conclude in Sec. IIXI 



II. ANDREEV REFLECTION 

The quantum mechanical description of Andreev re- 
flection starts from a pair of Schrodinger equations for 
electron and hole wave functions u(r) and u(r), coupled 
by the pair potential A(r). These socalled Bogoliubov- 
De Gennes (BdG) equations llj] take the form 



Hbg 



= E 



"^BG = 



H A(r) 
A*(r) ~H* 



(1) 
(2) 



The Hamiltonian H = {p + eA^ /2m + V - Ep is the 
single-electron Hamiltonian in the field of a vector po- 
tential A(r) and electrostatic potential V{y). The exci- 
tation energy E is measured relative to the Fermi energy 
Ep- If (w, v) is an eigenfunction with eigenvalue E, then 
{—v*,u*) is also an eigenfunction, with eigenvalue —E. 
The complete set of eigenvalues thus lies symmetrically 
around zero. The quasiparticle excitation spectrum con- 
sists of all positive E. 

In a uniform system with A(r) = A, A(r) = 0, V{y) = 
0, the solution of the BdG equations is 



E = [{h^k'^/2m-EFf + 1^^] 



1/2 



(3) 



u(r) == {2E)-^''^ {E + h^k^ I2m~ Epf^ e^^'^ (4) 



vi/2 ,k.> 



v{r) ^ {2E)-^'^E-h^ky2m + EF) e 



(5) 



The excitation spectrum is continuous, with excitation 
gap A. The eigenfunctions (u, v) are plane waves charac- 
terized by a wavevector k. The coefficients of the plane 
waves are the two coherence factors of the BCS (Bardeen- 
Cooper-Schrieffer) theory. 

At an interface between a normal metal and a su- 
perconductor the pairing interaction drops to zero over 
atomic distances at the normal side. (We assume non- 
interacting electrons in the normal region.) Therefore, 
A(r) = in the normal region. At the superconducting 
side of the NS interface, A(r) recovers its bulk value A 
only at some distance from the interface. This suppres- 
sion of A(r) is neglected in the step- function model 



A if res', 
if r e TV. 



(6) 



The step-function pair potential is also referred to in the 
literature as a "rigid boundary condition" jl^l ■ It greatly 
simplifies the analysis of the problem without changing 
the results in any qualitative way. 

Since we will only be considering a single superconduc- 
tor, the phase of the superconducting order parameter is 



irrelevant and we may take A real. (See Ref. H^ for a 
tutorial introduction to mesoscopic Josephson junctions, 
such as a quantum dot connected to two superconduc- 
tors.) 



III. MINIGAP IN NS JUNCTIONS 

The presence of a normal metal interface changes the 
excitation spectrum (proximity effect). The continuous 
spectrum above the bulk gap A differs from the BCS 
form (0J and in addition there may appear discrete en- 
ergy levels i?„ < A. 

The wave function of the lowest level contains elec- 
tron and hole components mq, vq of equal magnitude, 
mixed by Andreev reflection. The mean time TdwoU be- 
tween Andreev reflections (corresponding to the mean 
life time of an electron or hole excitation) sets the scale 
En = i?gap — ?i/Tdwcii for the energy of this lowest level 
This "minigap" is smaller than the bulk gap by a 



Eq = 



factor ^o/wi="''dwcii, with fo = Tiwf/A the superconducting 
coherence length and vp the Fermi velocity. The energy 
^/Tdweii is called the Thouless energy Et, because of the 
role it plays in Thouless's theory of localization Q. 

The simplest NS junction, which can be analyzed ex- 
actly jl5| | , consists of an impurity- free normal metal layer 
(thickness d) on top of a bulk superconductor. Because 
of translational invariance parallel to the NS interface, 
the parallel component pu of the momentum is a good 
quantum number. The lowest excitation energy 



EoiPw) 



2T{p« 



TiPii) = 



2md 



iP% 



■p2)l/2 



(7) 



is the reciprocal of the time T{p\\) between two subse- 
quent Andreev reflections. This time diverges when pn 
approaches the Fermi momentum pp = Tikp — \/2mEp^ 
so Eq can come microscopically close to zero. The lower 
limit Eq ^ h /md'^ is set by the quantization of the mo- 
mentum perpendicular to the layer. 

Impurities in the normal metal layer (with mean free 
path /) prevent the time between Andreev reflections to 
grow much lar ger than Ti j wsII — Ta.ax.(l / v p , d^ j v pi) . The 
excitation gap |T6lll7L ITsj 



E, 



h/Td„cii ^ (hvpl/d^) min(l, d^/f) (8) 



is now a factor kplinm{l,d^/P) larger than in the ab- 
sence of impurities. A precise calculation using disorder- 
averaged Green functions (reviewed in Ref. |ljj ) give s the 
curve shown in Fig. O The two asymptotes are [ig 



E — 



OAShvp/l, if d//< 1, 
0.78hD/d'^, if d/l:$> 1, 



(9) 



with D = vpl/S the diffusion constant in the normal 
metal. 

The minigap in a ballistic quantum dot (Andreev bil- 
liard) differs from that in a disordered NS junction in two 
qualitative ways: 
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FIG. 3: Excitation gap -Egap of a disordered NS junction, as 
a function of the ratio of the thickness d of the normal metal 
layer and the mean free path I. The curve in the bottom 
panel is calculated from the disorder- averaged Green function 
(for fo «C d,l). The top panel illustrates the geometry. The 
normal metal layer has a specularly reflecting upper surface 
and an ideally transmitting lower surface. Adapted from Ref. 
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1. The opening of an excitation gap depends on the 
shape of the boundary, rather than on the degree 
of disorder [23| . A chaotic billiard has a gap at the 
Thouless energy Et — h/rdwcii, like a disordered 
NS junction. An integrable billiard has a linearly 
vanishing density of states, like a ballistic NS junc- 
tion. 

2. In a chaotic billiard a new time scale appears, the 
Ehrenfest time Tp, which competes with Tdwcii hi 
setting the scale for the excitation gap |2j|. While 
■'■dwell is a classical ?i-independent time scale, Tp oc 
I ln?i| has a quantum mechanical origin. 

Because one can not perform a disorder average in 
Andreev billiards, the Green function formulation is less 
useful than in disordered NS junctions. Instead, we will 
make extensive use of the scattering matrix formulation, 
explained in the next section. 



IV. SCATTERING FORMULATION 

In the step-function model © the excitation spec- 
trum of the coupled electron-hole quasiparticles can be 
expressed entirely in terms of the scattering matrix of 
normal electrons [2^. 

The scattering geometry is illustrated in Fig. ^ It 
consists of a finite normal-metal region N adjacent to a 
semi- infinite superconducting region S. The metal region 
represents the Andreev billiard. To obtain a well-defined 
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FIG. 4: Normal metal (N) containing an Andreev billiard, 
coupled to a superconductor (S) by an ideal lead. The 
dashed line represents the NS interface. Scattering states 
c'" = {cf,c^) and c°"* = {c~,c^) are indicated schemati- 
cally. 



scattering problem we insert an ideal (impurity-free) nor- 
mal lead between N and S. We assume that the only 
scattering in the superconductor consists of Andreev re- 
flection at the NS interface (no disorder in S). The su- 
perconductor may then also be represented by an ideal 
lead. We choose a coordinate system so that the normal 
and superconducting leads lie along the x-axis, with the 
interface at a; = 0. 

We first construct a basis for the scattering matrix. In 
the normal lead N the eigenfunctions of the BdG equation 
(|T|) can be written in the form 



*«,c(N) = ( J ) ^ $„(y, z) exp(±ifc^x), (10a) 
0\ 1 



*„,h(N) 






$„(2/,z)exp(±ifcJja;),(10b) 



where the wavenumbers kf^ and fc„ are given by 



UG.h 



-iEF-E^a+(T'''''E) 



c,hJ^^l/2 



(11) 



and we have defined a° = 1, a^ = —1. The labels e 
and h indicate the electron or hole character of the wave 
function. The index n labels the modes, $„(2/, z) is the 
transverse wave function of the n-th mode, and £'„ its 
threshold energy: 



[ipl + pl)/2m + Viy, z)]^niy, z) = £;„$„(2/, z). (12) 



The eigenfunction <!>„ is normalized to unity, 
/d2//dz|$„|2 = l. 



In the superconducting lead S the eigenfunctions are 



*,tc(S) - 
*th(S) = 
We have defined 






1 



X $„(?/, z)exp(±iq,'; a; 

e''''V2 \ 1 
e-'^'/s 



(£;2/a2 _ 1) 



-1/4 



(13a) 



X $„(y,z)exp(±i(jj;x 



(£;Va2 - 1) 



-1/4 



(13b) 



[Ef - En + fj'='^(£;2 - A2)1/2]1/2^ (^4) 



v 



■^ = CT"''^arccos(£;/A). 



(15) 



The wave functions pO|l and p3(l have been normal- 
ized to carry the same amount of quasiparticle current, 
because we want to use them as the basis for a unitary 
scattering matrix. The direction of the velocity is the 
same as the wave vector for the electron and opposite for 
the hole. 

A wave incident on the Andreev billiard is described 
in the basis UlUfl by a vector of coefficients 



(4>Ch), 



(16) 



as shown schematically in Fig. ^ (The mode index n has 
been suppressed for simplicity of notation.) The reflected 
wave has vector of coefficients 



(Ce-,C+) 



(17) 



The scattering matrix iSn of the normal region relates 
these two vectors, c^"* = ^ncJ^jV Because the normal 
region does not couple electrons and holes, this matrix 
has the block-diagonal form 



S^{E) 



( S{E) 

(^ s{-Ey 



(18) 



Here S{E) is the unitary scattering matrix associated 
with the single-electron Hamiltonian H . It is an TV x A^ 
matrix, with N{E) the number of propagating modes at 
energy E. The dimension of S'n(^) is N{E) + N{-E). 

For energies < -E < A there are no propagating 
modes in the superconducting lead S. Restricting our- 
selves to that energy range, we can define a scattering 
matrix Sa for Andreev reflection at the NS interface by 
c"^ = Sac°^^. The elements of Sa are obtained by match- 
ing the wave function (|10l) at a; = to the decaying wave 
function (|13|) . Since A <^ Ep one may ignore normal 
reflections at the NS interface and neglect the difference 
between N{E) and Ni—E). This is known as the An- 
dreev approximation []|. The result is 



Sa{E) = 



a{E) 



a{E) 
a{E) 



-i arccos(£'/A) 



E 

A 



i\ I 



_E2 
A2- 



(19) 
(20) 



Andreev reflection transforms an electron mode into a 
hole mode, without change of mode index. The transfor- 
mation is accompanied by a phase shift — arccos(-B/A) 
due to the penetration of the wave function into the su- 
perconductor. 

We are now ready to relate the excitation spectrum 
of the Andreev billiard to the scattering matrix of the 
normal region. We restrict ourselves to the discrete 
spectrum (see Ref. |23 for the continuous spectrum). 
The condition Qn = ^A^NCin for a bound state implies 
Det (1 - SaSn) = 0. Using Eqs. (0, (0, and the iden- 
tity 

Det f ^ ^ W Det {ad ~ aca-^b) (21) 

one obtains the equation |22| 

Det [l - a{EfS{E)S{~-E)*] = 0. (22) 

The roots Ep of this determinantal equation constitute 
the discrete spectrum of the Andreev billiard. 



V. STROBOSCOPIC MODEL 

Although the phase space of the Andreev billiard is 
four-dimensional, like for any billiard it can be reduced 
to two dimensions on a Poincare surface of section |3, Q ■ 
This amounts to a stroboscopic description of the clas- 
sical dynamics, because the position and momentum are 
only recorded when the particle crosses the surface of sec- 
tion. Quantum mechanically, the stroboscopic evolution 
of the wave function is described by a compact unitary 
map rather than by a noncompact Hcrmitian operator 
[2a. |2J| . What one loses by the stroboscopic description 
is information on time scales below the time of flight 
across the billiard. What one gains is an enormous in- 
crease in computational efficiency. 

A stroboscopic model of an Andreev billiard was con- 
structed by Jacquod et al. l25|, building on an existing 
model for open normal billiards called the open kicked 
rotator J2y|. The Andreev kicked rotator possesses the 
same phenomenology as the Andreev billiard, but is much 
more tractable numerically.^ In this subsection we dis- 
cuss how it is formulated. Some results obtained by this 
numerical method will be compared in subsequent sec- 
tions with results obtained by analytical means. 

A compact unitary map is represented in quantum 
mechanics by the Floquet operator F, which gives the 
stroboscopic time evolution u{pTo) ~ F''w(0) of an ini- 
tial wave function m(0). (We set the stroboscopic period 
To = 1 in most equations.) The unitary M x M ma- 
trix F has eigenvalues exp(— iem), with the quasi-energies 



^ The largest simulation to date of a two-dimensional Andreev 
billiard has N = 30, while for the Andreev kicked rotator N = 
10^ is within reach, cf. Fig. |^ 



Sm S (— TT, tt) (measured in units of h/ro). This describes 
the electron excitations above the Fermi level. Hole exci- 
tations below the Fermi level have Floquet operator F* 
and wave function v{p) = (F*)Pw(0). The mean level 
spacing of electrons and holes separately is (5 = 2tt/M. 

An electron is converted into a hole by Andreev re- 
flection at the NS interface, with phase shift —i for 
£ ^ ToA/h [cf. Eq. l(^ ]. In the stroboscopic descrip- 
tion one assumes that Andreev reflection occurs only at 
times which are multiples of tq. The N x M matrix P 
projects onto the NS interface. Its elements are P„m = 1 
if m = n e {rii, n2, . . . n^} and Pnm = otherwise. The 
dwell time of a quasiparticle excitation in the normal 
metal is TdwcU — M /N , equal to the mean time between 
Andreev reflections. 

Putting all this together one constructs the quantum 
Andreev map from the matrix product 



T = V 



F 

F* 



V = 



1 - P^P -iP^P 
-iP^P 1 - pTp 



(23) 

(The superscript "T" indicates the transpose of a ma- 
trix.) The particle-hole wave function \E' = {u,v) evolves 
in time as ^(p) = J^^^(O). The Floquet operator can be 
symmetrized (without changing its eigenvalues) by the 
unitary transformation JF -^ 'p-^/'^j^'p'^/^^ with 



7?i/2 



1 - (1 - iV2)P^P -i^V2P^P 



i^V2P^P l-(l-i^/2)FTp 



(24) 
The quantization condition det{J-' — e *'^) = can be 
written equivalently as [25J 

Det [1 + S'(e)S'(-e)*] = 0, (25) 

in terms of the N x N scattering matrix [2^ |23 

(26) 



S{e) = Pie-'" - F(l - P^P)]-^FP'^ 



Eq. (|25|l for the Andreev map has the same form as Eq. 
(|22|l for the Andreev billiard (with a -^ —i). In par- 
ticular, both equations have roots that lie symmetrically 
around zero. 

A specific realization of the Andreev map is the An- 
dreev kicked rotator. (See Ref. [23| for a different real- 
ization, based on the kicked Harper model.) The normal 
kicked rotator has Floquet operator [23 



'"P I *4io 9^ ' '"P 



KIo , 
-I— — cos( 

TlTr, 



X exp 






(27) 



It describes a particle that moves freely along the unit 
circle (cos 0, sin 6*) with moment of inertia Iq for half a 
period tq, is then kicked with a strength Kcos9, and 
proceeds freely for another half period. Upon increas- 
ing K the classical dynamics varies from fully integrable 



{K = 0) to fully chaotic [K > 7, with Lyapunov ex- 
ponent a « ln(i^/2)]. For K < 7 stable and unstable 
motion coexist (mixed phase space). If needed, a mag- 
netic field can be introduced into the model as described 
inRef. jS^- 

The transition from classical to quantum behavior 
is governed by the effective Planck constant hcff = 
fiTo/2TrIo. For 1/hcS = M an even integer, F can be 
represented by an Af x M unitary symmetric matrix. 
The angular coordinate and momentum eigenvalues are 
9m = 2T:m/M and Pm — fim, with m = 1,2,...M, so 
phase space has the topology of a torus. The NS in- 
terface is an annulus around the torus, either in the 9- 
direction or in the p-direction. (The two configurations 
give equivalent results.) The construction H23(l produces 
a 2M X 2M Floquet operator JF, which can be diago- 
nalized efHciently in 0{M'^ In Af ) operations [rather than 
0{M^)] by combining the Lanczos technique with the 
fast-Fourier-transform algorithm jsj . 



VI. RANDOM-MATRIX THEORY 

An ensemble of isolated chaotic billiards, constructed 
by varying the shape at constant area, corresponds to 
an ensemble of Hamiltonians H with a particular dis- 
tribution function P{H). It is convenient to think of 
the Hamiltonian as a random M x M Hermitian ma- 
trix, eventually sending M to infinity. The basic pos- 
tulate of random-matrix theory (RMT) 9] is that the 
distribution is invariant under the unitary transforma- 
tion H — > UHU\ with U an arbitrary unitary matrix. 
This implies a distribution of the form 



P{H) cxexp[-Try(iy)]. 



(28) 



If V{H) oc iJ^, the ensemble is called Gaussian. This 
choice simplifies some of the calculations but is not es- 
sential, because the spectral correlations become largely 
independent of V in the limit Af — > oo. More generally, 
the ensemble of the form H28|) is called the Wigner-Dyson 
ensemble, after the founding fathers of RMT. 

By computing the Jacobian from the space of matrix 
elements to the space of eigenvalues En (n = 1, 2, . . . M), 
one obtains the eigenvalue probability distribution Q 



P{{E,,})^\{\E.,-E,f\{^ 



-V{Ek) 



(29) 



Kj 



The symmetry index (3 counts the number of degrees of 
freedom in the matrix elements. These are real (/3 = 1) 
in the presence of time-reversal symmetry or complex 
(/3 = 2) in its absence. (A third possibility, /3 = 4, ap- 
plies to time-reversally symmetric systems with strong 
spin-orbit scattering, which we will not consider here.) 
Since the unitary transformation H — > UHW requires 
an orthogonal U to keep a real Hamiltonian, one speaks 
of the Gaussian orthogonal ensemble (GOE) when /3 = 1. 




FIG. 5: A quantum dot (N) connected to a superconductor 
(S). The voltages on the gates Vi and V2 change the shape 
of the dot. Different values of the applied voltages create 
different samples within the same ensemble. From Ref. |3j|. 



The name Gaussian unitary ensemble (GUE) refers to 
/3 = 2. ^ 

There is overwhelming numerical evidence that chaotic 
billiards are well described by the Wigner-Dyson ensem- 
ble [3. (This is known as the Bohigas-Giannoni-Schmit 
conjecture |33|.) A complete theoretical justification is 
still lacking, but much progress has been made in that 
direction |33|- In this section we will take Eq. H28() for 
the ensemble of isolated billiards as our starting point 
and deduce what properties it implies for the ensemble 
of Andreev billiards. 

The isolated billiard becomes an Andreev billiard when 
it is connected by a point contact to a superconductor, 
cf. Fig. [5| In the isolated billiard RMT breaks down on 
energy scales greater than Ti/tc-c^, with the ergodic time 
Torg — A^/"^ /vp set by the time of flight across the bil- 
liard (of area A, at Fermi velocity vp). On larger energy 
scales, hence on shorter time scales, non-chaotic dynam- 
ics appears which is beyond RMT. The superconductor 
affects the billiard in an energy range around the Fermi 
level that is set by the Thouless energy Et — ?i/Tdwcn- 
(We assume that Et is less than the gap A in the bulk 
superconductor.) In this context the dwell time Tdwcii is 
the mean time between Andreev reflections (being the 
life time of an electron or hole quasiparticle) . The condi- 
tion Torg <C Tdwoii of weak coupling is therefore sufflcient 
to be able to apply RMT to the entire relevant energy 
range. 



A. Effective Hamiltonian 

The excitation energies Ep of the Andreev billiard in 
the discrete part of the spectrum are the solutions of the 
determinantal equation (|22l) , given in terms of the scat- 
tering matrix S{E) in the normal state (i.e. when the su- 
perconductor is replaced by a normal metal) . This equa- 
tion can alternatively be written in terms of the Hamil- 
tonian H of the isolated billiard and the M x N coupling 
matrix W that describes the A^-mode point contact. The 



relation between S and if , VF is [j, nU^ 

S{E) = l-27:iW^{E-H + iTTWW'^)-'^W. (30) 
The N X N matrix W^W has eigenvalues w„ given by 

wn = -^ (2 - r„ - 2v/i - r„) , (31) 

where d is the mean level spacing in the isolated billiard 
and r„ G [0, 1] is the transmission probability of mode 
n = 1,2, ... N in the point contact. For a ballistic con- 
tact, r„ = 1, while r„ ^ 1 for a tunneling contact. Both 
the number of modes N and the level spacing S refer to 
a single spin direction. 

Substituting Eq. (|^ into Eq. (|^ . one arrives at an 
alternative determinantal equation for the discrete spec- 
trum 35] : 

Det [E-n + W{E)] = 0, (32) 

—H* I ' ' 

EWW'^ AWW^ 



The density of states follows from 



(34) 



P{E) 



1 



lmTi{l+dW/dE){E+iO+-n+Wy\ (35) 



In the relevant energy range E < Et <C A the ma- 
trix yV{E) becomes energy independent. The excitation 
energies can then be obtained as the eigenvalues of the 
effective Hamiltonian ISal 



'HeS — 



H -ttWW'^ 

-ttWW^ -H* 



(36) 



The effective Hamiltonian Hcff should not be confused 
with the Bogoliubov-de Gennes Hamiltonian TYbg i which 
contains the superconducting order parameter in the off- 
diagonal blocks [cf. Eq. 0]. The Hamiltonian Hbg de- 
termines the entire excitation spectrum (both the dis- 
crete part below A and the continuous part above A), 
while the effective Hamiltonian Tieff determines only the 
low-lying excitations Ep <^ A. 

The Hermitian matrix Tieff (like TYbg ) is antisymmet- 
ric under the combined operation of charge conjugation 
(C) and time inversion (T) [37] : 



TicS = —^yT^cS^yi ^y 







(37) 



(An M X M unit matrix in each of the four blocks of ay is 
implicit.) The CT-antisymmetry ensures that the eigen- 
values lie symmetrically around E = 0. Only the positive 
eigenvalues are retained in the excitation spectrum, but 
the presence of the negative eigenvalues is felt as a level 
repulsion near E = 0. 



B. Excitation gap 

In zero magnetic field the suppression of the density of 
states p{E) around E = extends over an energy range 
Et that may contain many level spacings S of the iso- 
lated billiard. The ratio g ~ Et/S is the conductance of 
the point contact in units of the conductance quantum 
e^//i. For g ^ 1 the excitation gap fi'gap — 5^ is a meso- 
scopic quantity, because it is intermediate between the 
microscopic energy scale S and the macroscopic energy 
scale A. One can use perturbation theory in the small 
parameter 1/g to calculate p{E). The analysis presented 
here follows the RMT of Melsen et al. i20(]. An alter- 
native derivation [Sg , using the disorder-averaged Green 
function, is discussed in the next sub-section. 

In the presence of time-reversal symmetry the Hamil- 
tonian H of the isolated billiard is a real symmetric ma- 
trix. The appropriate RMT ensemble is the GOE, with 
distribution ^] 



P{H) ex exp 



AMS^ 



TtH^ 



(38) 



The ensemble average (• • • ) is an average over H in the 
GOE at fixed coupling matrix W. Because of the block 
structure of Ticff , the ensemble averaged Green function 
g{E) = {{E - Hcfty^) consists of four M x M blocks 
Gil, G12, G21, G22- By taking the trace of each block 
separately, one arrives at a 2 x 2 matrix Green function 



G 



Gu G12 
G21 G22 



TrtJii TtGi2 
TTG21 TtG22 



(39) 



(The factor S/t: is inserted for later convenience.) 

The average over the distribution 138|) can be done 
diagrammatically 39, 40]. To leading order in 1/AI and 
ior E ^ S only simple (planar) diagrams need to be 
considered. Resummation of these diagrams leads to the 
selfconsistency equation [20I 1351 ] 



G = [E+W -{MS /Tr)a^Ga^]-\ a. 



1 
-1 



(40) 



This is a matrix-generalization of Pastur's equation in 
the RMT of normal systems 41]. 

The matrices in Eq. l|in|l have four M x M blocks. By 
taking the trace of each block one obtains an equation 
for a 4 X 4 matrix, 



G 



M 

M ^ 

771—1 



TXEJMS - Gil Wm + G12 

wm + G21 t:E/MS - G22 



TT^w^/MS if TO= l,2,...iV, 
if TO = A^ + 1, . . . M. 



(4^) 
(42) 



Since G22 — Gu and G21 — G12 there are two unknown 
functions to determine. For M ^ N these satisfy 

G^2 = 1 + G^i, 



—^Gi2 - Gii2j( 



(43a) 
-G12 + 1 - 2/r„)-\ (43b) 
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FIG. 6: Ensemble averaged density of states of a chaotic bil- 
liard coupled by a point contact to a superconductor, for 
several values of the transmission probability through the 
point contact. The energy is in units of the Thouless energy 
Et ~ NT5/4:TT. The solid curves are computed from Eqs. 
(I43I I and 14411 . for mode-independent transmission probabili- 
ties r = 1, 0.5, 0.25, 0.1. The dashed curve is the asymptotic 
result (EU for r < 1. Adapted from Ref. |20|. (The definition 
of 5 used in that paper differs from the one used here by a 
factor of two.) 



A rather simple closed-form expression for {p{E)) ex- 
ists in two limiting cases |23- In the case F = 1 of a 
ballistic point contact one has 

{piE)) = ^^[Q+{E/Et) - Q-iE/Er)], (47) 



3E6 



-ix) 



8 - 362:^ ± 3x^3x4 + 132x2 - 48 
0.30 ;i 



E > Eg^p = 2-i'>/'Et = OMEt 



TdwcU 



1/3 

,(48) 

0.048 7V(5, 
(49) 



where 7 = ■^{yE — 1) is the golden number. In this case 
the parameter Agap in Eq. H46|l is given by 

Agap = [(5 - 2 V5)(52£;gap/87r2] ^/^ ^ 0.O68 N^'^5. (50) 
In the opposite tunneling limit F ^ 1 one finds 



{p{E)) = ^{E' - El)-^'\ E > i^gap = Et. 



(51) 



In this limit the density of states of the Andreev billiard 
has the same form as in the BCS theory for a bulk super- 
conductor |42| , with a reduced value of the gap ( "mini- 
gap"). The inverse square-root singularity at the gap is 
cut off for any finite F, cf. Fig. 



where we have used the relation (|31|l between the pa- 
rameters Wn and the transmission probabilities F„. Eq. 
H43|l has multiple solutions. The physical solution satis- 
fies l[mE^oo{p{E)) = 2/(5, when substituted into 



(p(^))=-(2/<5)ImGii(^). 



(44) 



In Fig. El we plot the density of states in the mode- 
independent case F„ = F, for several vaues of F. It van- 
ishes as a square root near the excitation gap. The value 
of Egs,p can be determined directly by solving Eq. l(^ 
jointly with dE/dGu = 0. The result is 



(1 - ky 



3k^ - 20k^ 



16 



[i-ky 



3fc2 + 8 
(l-fc)2' 



X = Eg^p/Er, k = l- 2/F, Et ^ NTS/An. 



= 1, 
(45) 



For later use we parametrize the square-root dependence 
near the gap as 



{piE)) ^ - 

IT 



1 E-E, 



:ap 



A3 
gap 



E~^ E, 



gap- 



(46) 



When E ^ ^gap the density of states approaches the 
value 2/S from above, twice the value in the isolated bil- 
liard. The doubling of the density of states occurs be- 
cause electron and hole excitations are combined in the 
excitation spectrum of the Andreev billiard, while in an 
isolated billiard electron and hole excitations are consid- 
ered separately. 



C. Effect of impurity scattering 

Impurity scattering in a chaotic Andreev billiard re- 
duces the magnitude of the excitation gap by increasing 
the mean time Td^oii between Andreev reflections. This 
effect was calculated by Vavilov and Larkin [Sg using the 
method of impurity-averaged Green functions [l9j. The 
minigap in a disordered quantum dot is qualitatively sim- 
ilar to that in a disordered NS junction, cf. Sec. lIIII The 
main parameter is the ratio of the mean free path I and 
the width of the contact W. (We assume that there is no 
barrier in the point contact, otherwise the tunnel proba- 
bility F would enter as well.) 

For I ^ W the mean dwell time saturates at the bal- 
listic value 



T'dwcll 



2Trn 

'ns 



■kA 



if 1~>W. 



(52) 



In the opposite limite / <C W the mean dwell time is 
determined by the two-dimensional diffusion equation. 
Up to a geometry-dependent coefficient c of order unity, 
one has 



T-dwcii - c—^ ln(A/VF2)^ if ; ^ ^_ 



(53) 



The density of states in the two limits is shown in Fig. 
13 There is little difference, once the energy is scaled 
by TdwcU- For / ^ W the excitation gap is given by 



the RMT result E, 



ap 



0.300 n/Td 



well 



I <S^ W Vavilov and Larkin find E,y 



cf. Eq. (gni 
: 0.331 ?i/ TdwcU 



For 




dwell 



FIG. 7: The dashed curve is the F = 1 result of Fig. |S1 cor- 
responding to a quantum dot with weak impurity scatter- 
ing (mean free path I much larger than the width W of the 
point contact). The solid curve is the corresponding result for 
strong impurity scattering [l <C W). The line shape is almost 
the same, but the energy scale is different (given by Eqs. H52|l 
and 15311 ■ respectively). Adapted from Ref. |3a|. 



D. Magnetic field dependence 

A magnetic field B, perpendicular to the billiard, 
breaks time-reversal symmetry, thereby suppressing the 
excitation gap. A perturbative treatment remains possi- 
ble as long as £'gap(i3) remains large compared to 5 J43l |. 

The appropriate RMT ensemble for the isolated bil- 
liard is described by the Pandey-Mehta distribution 



P{H) ex cxp 



71^(1 +b^) 






(54) 



The parameter h E [0, 1] measures the strength of the 
time-reversal symmetry breaking. The invariance of 
P{H) under unitary transformations is broken if 6 ^^ 0, T 
The relation between b and the magnetic flux $ through 
the billiard is Q 



Mb"^ = c{^e/h)- 



hvp 



(55) 



with c a numerical coefficient that depends only on the 
shape of the billiard. Time-reversal symmetry is effec- 
tively broken when Mb^ ~ g, which occurs for $ ~ 
(/i/e)-\/Tcrg/rdweii "C h/e. The effect of such weak mag- 
netic fields on the bulk superconductor can be ignored. 

The selfconsistency equation for the Green function is 
the same as Eq. (|41|) . with one difference: On the right- 
hand-side the terms G12 and G21 are multiplied by the 
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FIG. 8: Magnetic field dependence of the density of states for 
the case of a ballistic point contact (r„ = 1), computed from 
Eqs. H43a^ . H44^ . and 1)56^ . The microscopic gap of order 5 
which persists when $ > <I>c is not resolved in this calculation. 
Adapted from Ref. |4^ . 



factor (l-fe2)/(i + 62). In the limit M ^ cx) , b ^ 0, Mb"^ 
finite, the first equation H43a|l still holds, but the second 
equation Ij43b|) is replaced by 



{2ttE/5 - 4M62G'ii)G'i2 = Gn 

AT 



(56) 



The resulting magnetic field dependence of the average 
density of states is plotted in Fig. (SJ for the case r„ = 1 
of a ballistic point contact. The gap closes when Mb"^ = 
A^r/S. The corresponding critical flux $c follows from 
Eq. |5^. 



E. Broken time-reversal symmetry 

A microscopic suppression of the density of states 
around E = Q, on an energy scale of the order of the 
level spacing, persists even if time-reversal symmetry is 
fully broken. The suppression is a consequence of the 
level repulsion between the lowest excitation energy Ei 
and its mirror image —Ei, which itself follows from the 
CT-antisymmetry (|37|) of the Hamiltonian. Because of 
this mirror symmetry, the effective Hamiltonian Tieft of 
the Andreev billiard can be factorized as 



'HcS — U 



£ 
-£ 



U\ 



(57) 



with U a 2M x 2M unitary matrix and £ — 
diag{Ei, E2, ■■■ Em) a diagonal matrix containing the 
positive excitation energies. 

Altland and Zirnbauer [37| have surmised that an en- 
semble of Andreev billiards in a strong magnetic field 



10 



would have a distribution of Hamiltonians of the Wigner- 
Dyson form (|28|l . constrained by Eq. H57I) . This con- 
straint changes the Jacobian from the space of matrix 
elements to the space of eigenvalues, so that the eigen- 
value probability distribution is changed from the form 
(|^ (with /3 = 2) into 



Pi{E.n})^l[iEf-Effl[Ele 



-V{Ek)-V{-Ek) 



(58) 



i<j 



The distribution ^^ with V{E) ex E'^ is related to 
the Laguerre unitary ensemble (LUE) of RMT 9] by a 
change of variables. The average density of states van- 
ishes quadratically near zero energy [45J, 



, , ,, ^ / iimUTrE/S)\ 



(59) 



All of this is qualitatively different from the "folded 
GUE" that one would obtain by simply combining two 
independent GUE's of electrons and holes |4fil |. 

A derivation of Altland and Zirnbauer's surmise has 
been given by Frahm et al. [Sg, who showed that the 
LUE for the effective Hamiltonian Ticff of the Andreev 
billiard follows from the GUE for the Hamiltonian H of 
the isolated billiard, provided that the coupling to the 
superconductor is sufficiently strong. To compute the 
spectral statistics on the scale of the level spacing, a non- 
perturbative technique is needed. This is provided by 
the supersymmetric method jig. (See Ref. ^SJ for an 
alternative approach using quantum graphs.) 

The resulting average density of states is [3g 



{p{E)) = - 



2 sm{2TTE/5) 
'5 ^ 



ds 







/2ttE / 4J 

X cos — — W IH 

V '' V 9A 



9A 



N 

E 



2r2 



1(2 



(60) 
(61) 



The parameter gA is the Andreev conductance of the 
point contact that couples the billiard to the supercon- 
ductor 1431. The Andreev conductance can be much 



N 



smaller than the normal-state conductance g = X]n=i ^n- 
(Both conductances are in units of 2e^//i.) In the tunnel- 
ing hmit r„ = r < 1 one has g = NT while ffA = 5 NT"^. 
Eq. © describes the crossover from the GUE result 
p{E) = 2/S for gA < 1 to the LUE result ^ for 
(7A ^ 1. The opening of the gap as the coupling to 
the superconductor is increased is plotted in Fig. El The 
CT-antisymmetry becomes effective at an energy E for 
^A <] E/S. For small energies E <^ S min{y/gX, 1) the 
density of states vanishes quadratically, regardless of how 
weak the coupling is. 
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FIG. 9: Density of states of the Andreev billiard in a strong 
magnetic field for three different values of the Andreev con- 
ductance of the point contact: gA ~ 0.4, 4, 40. The solid 
curves are calculated from Eq. HtJUl . The dashed line is the 
LUE result 1591 . corresponding to the limit (;a ^ 00. The 
dotted line is the GUE limit qa — > 0. Adapted from Ref. ^g]. 



F. Mesoscopic fluctuations of the gap 

The smallest excitation energy Ei in the Andreev bil- 
liard fluctuates from one member of the ensemble to the 
other. Vavilov et al. [sj have surmised that the dis- 
tribution of these fluctuations is identical upon rescal- 
ing to the known distribution 50] of the lowest eigen- 
value in the Gaussian ensembles of RMT. This surmise 
was proven using the supersymmetry technique by Os- 
trovsky, Skvortsov, and Feigelman I5l| and by Lamacraft 
and Simons 52] . Rescaling amounts to a change of vari- 
ables from El to X — {Ei — £^gap)/Agap, where Sgap and 
Agap parameterize the square-root dependence (|3^ of 
the mean density of states near the gap in perturbation 
theory. The gap fluctuations are a mesoscopic, rather 
than a microscopic effect, because the typical magnitude 
Agap ~ i?gap(5^/^ of the fluctuations is > (5 for Eg^p > S. 
Still, the fluctuations are small on the scale of the gap 
itself. 

Following Ref. [SJI , in zero magnetic field the gap dis- 
tribution is obtained by rescaling the GOE result of Tracy 
and Widom [Hfl , 



^(^l) = ^^1 [(El ~ i?gap)/Agap] , 

Fi{x) = cxp ( -i f [q{x') + {x- x')q^{x')]dx' 



(62) 



The function q{x) is the solution of 

q"{x) = -xq{x) + 2q^{x), 



(63) 



(64) 



with asymptotic behavior q{x) — > Ai(— x) as x — > — cx) 
[Ai(x) being the Airy function]. For small x there is a 
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Energy scale 



Flux scale 



FIG. 10: Probability distribution of the rescaled excitation 
gap X — (El — iJgap)/Agap, in the presence [/3 — 1, Eq. (|S5J] 
and absence [f3 = 2, Eq. lISZH of time-reversal symmetry. 
Adapted from Ref. ^. 



tail of the form 



Pix) 



4V^|x|i/4 



exp 



(-11-1 



3/2 



a;< -1. 



(65) 



The distribution (|62|l is shown in Fig. ^] (solid curve) . 
The mean and standard deviation are 



(El) = ^gap + 1.21 Agap, {{El - {Ei))y/^ = 1.27 A, 



gap- 



Because the mesoscopic fluctuations in the gap occur 
on a much smaller energy scale Agap than ii'gap, there 
exists a range of magnetic fields that break time-reversal 
symmetry of the gap fluctuations without significantly 
reducing E^gap |3J|. In this field range, specified in Table 
HI the distribution of the lowest excitation is given by the 
GUE result M 



P{El) = ^mEl - i?gap)/Agap], (67) 



F2 {x) exp - / {x- x')q^ {x')dx' 



(68) 



This curve is shown dashed in Fig. ^| The tail for small 
x is now given by 



P{x) 



1 



SttIiI 



exp 



(-II-I 



3/2 



, a;< -1. 



(69) 



The mesoscopic gap fluctuations induce a tail in the 
ensemble averaged density of states {p{E)) for E < E^^. 
In the same rescaled variable x the tail is given by [sj 

{p{x)) = ~xA\^{x) + [M'{x)Y 

/"OO 

+ \5p,xM{x)[l~ Ai{y)dy]. (70) 



Bulk statistics 

Edge statistics 

Gap size 



£.1/3^2/3 



{h/e 



1/2 



s^^^/h^ 



/2 



(/./e)ri/^5i/6^i/3/^i/2 

(Ve)TcVg'Sga/p/?l'''' 



TABLE I: Characteristic energy and magnetic flux scales 
for the spectral statistics in the bulk and at the edge of the 
spectrum and for the size of the gap. The (3 = 2 distribution 
l| W| l applies to the flux range {h/e)Ti^^S^/'^Elip/h^/^ < $ < 
(Ve)r^/g'£^/^/fti/2. 




FIG. 11: Ensemble averaged density of states (p) together 

with the probability distribution P of the excitation gap, as 

(66 j a function of the rescaled energy x = {E — _Bgap)/Agap. The 



dotted and dashed curves are the universal results II62|I and 
I70II of RMT in the presence of time-reversal symmetry (/3 
1). The solid curve is the mean density of states 14611 
perturbation theory. Adapted from Ref. 



Asymptotically, {p{x)) ex exp(— |/3|a;p/^) for x <C —1. 
The tail in {p{E)) is the same as the tail in P{E), as it 
should be, since both tails are due to the lowest eigen- 
value. In Fig. Illl we compare these two functions in zero 
magnetic field, together with the square-root density of 
states from perturbation theory. 

A numerical simulation of the stroboscopic model of 
Sec. IVl provides a test of these predictions I23. Results 
are shown in Fig. ^1 for the case (3=1 and deep in 
the chaotic regime (kicking strength K ^ 1). The agree- 
ment with RMT is very good — without any adjustable 
parameters. 



G. Coulomb blockade 

Coulomb interactions between electron and hole quasi- 
particles break the charge-conjugation invariance H37|) of 
the Hamiltonian. Since Andreev reflection changes the 
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FIG. 12; Main plot: Gap distribution for the Andreev kicked 
rotator with parameters M = 2n/5 — 8192, kicking strength 
K = 45, and M/N = rd„oii = 10 (o), 20 (•), 40 (+), and 
50 (x). There is no magnetic field. The solid line is the 
RMT prediction 16211 . Inset: Average density of states for 
the same system. The solid line is the RMT prediction (|49^ . 
(Deviations from perturbation theory are not visible on the 
scale of the inset.) Adapted from Ref. J2RJ |. 



charge on the billiard by 2e, this scattering process be- 
comes energetically unfavorable if the charging energy 
Ec exceeds the superconducting condensation energy 
(Josephson energy) Ej. For Ec ^ Ej one obtains the 
Coulomb blockade of the proximity effect studied by Os- 
trovsky, Skvortsov, and Feigelman [S^. 

The charging energy Ec = e^/2C is determined by the 
capacitance C of the billiard. The Josephson energy is 
determined by the change in free energy of the billiard 
resulting from the coupling to the superconductor, 



E, = - 



[p{E) - 2/(5] EdE. 



(71) 



The discrete spectrum E < i?gap contributes an amount 
of order E'l /5 to Ej. In the continuous spectrum 
E > i^gap the density of states p{E), calculated by RMT, 
decays ex l/E"^ to its asymptotic value 2/5. This leads to 
a logarithmic divergence of the Josephson energy ,35..54| . 
with a cutoff set by min(A, U/Tag): 



E,^ 



E. 



gap j^ I min(A, ?i/Torg) 



E, 



:ap 



(72) 



The suppression of the excitation gap with increasing 
En is plotted in Fig. [121 for the case F < 1, A < h/Tcrg 
[53J. The initial decay is a square root. 



1 - Aeff/Sg 



EcS 



EI,^H2A/Eg,p) 
and the final decay is exponential. 



1/2 



«1, (73) 




Aeff/A = 2exp(-2^c'5/^|ap)«l. (74) 



EcS/E^ap 

FIG. 13: Suppression due to Coulomb interactions of the gap 
Aeff in the density of states of an Andreev billiard coupled by 
a tunnel junction to a superconductor, relative to the nonin- 
teracting gap Sgap = NTS /in (with F < 1 < iVF). The plot 
is for the case A = e^E^^p <^ fi/ra-c^. The dashed lines are 
the asymptotes |I7^ and lltHl . Adapted from Ref. [H^. 



Here Aoff refers to the gap in the presence of Coulomb 
interactions and -Egap = NV5 /A.-K is the noninteracting 
value (EB- 

The gap Acff governs the thermodynamic properties of 
the Andreev billiard, most importantly the critical cur- 
rent. It is not, however, the relevant energy scale for 
transport properties. Injection of charge into the bil- 
liard via a separate tunnel contact measures the tunnel- 
ing density of states ptunnci, which differs in the presence 
of Coulomb interactions from the thermodynamic den- 
sity of states p considered so far. The gap Atunnci in 
Ptunnei crosscs ovcr from the proximity gap i?gap when 
Ec <C Ej to the Coulomb gap Ec when Ec > Ej, see 
Fig. El The single peak in ptunnoi at Atunnci splits into 
two peaks when Ec and Ej are of comparable magnitude 
[53 ■ This peak splitting happens because two states of 
charge -\-e and — e having the same charging energy are 
mixed by Andreev reflection into symmetric and anti- 
symmetric linear combinations with a slightly different 
energy. 



VII. QUASICLASSICAL THEORY 

It was noticed by Kosztin, Maslov, and Goldbart p 
that the classical dynamics at the Fermi energy in an 
Andreev billiard is integrable — even if the dynamics in 
the isolated billiard is chaotic. Andreev reflection sup- 
presses chaotic dynamics because it introduces a peri- 
odicity into the orbits: The trajectory of an electron is 
retraced by the Andreev reflected hole. At the Fermi en- 
ergy the hole is precisely the time reverse of the electron, 
so that the motion is strictly periodic. For finite excita- 
tion energy or in a non-zero magnetic field the electron 
and the hole follow slightly different trajectories, so the 
orbit does not quite close and drifts around in phase space 
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FIG. 14: Main plot: gap Atunnci in the tunneling density of 
states as a function of the charging energy (for A = e^Eg^p 
and A'^r = 40 tt). The initial decay (barely visible on the 
scale of the plot) follows Eq. 17311 and crosses over to an in- 
crease (Atunnci ^^ Ec)- Inset: tunneling density of states at 
EcS/Eg^p = 2.8 (corresponding to the dot in the main plot). 



Adapted from Ref. ,53.1 . 
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The near-periodicity of the orbits implies the existence 
of an adiabatic invariant. Quantization of this invariant 
leads to the quasiclassical theory of Silvestrov et al. [53 ■ 



A. Adiabatic quantization 




FIG. 15: Classical trajectory in an Andreev billiard. Particles 
in a two-dimensional electron gas are deflected by an elec- 
trostatic potential. (The dotted circles are equipotentials.) 
There is specular reflection at the boundaries with an insula- 
tor (thick solid lines) and Andreev reflection at the boundary 
with a superconductor (dashed line). The trajectory follows 
the motion between two Andreev reflections of an electron 
near the Fermi energy. The Andreev reflected hole retraces 
this trajectory in opposite direction. From Ref. |59| . 



Figs. ll5| and ll6l illustrate the nearly periodic motion in 
a particular Andreev billiard. Fig. 1151 shows a trajectory 
in real space while Fig. ^| is a section of phase space 
at the interface with the superconductor (y — 0). The 
tangential component Px of the electron momentum is 
plotted as a function of the coordinate x along the inter- 
face. Each point in this Poincare map corresponds to one 
collision of an electron with the interface. (The collisions 
of holes are not plotted.) The electron is retroreflected as 
a hole with the same p^- At the Fermi level {E = 0) the 
component Py is also the same, and so the hole retraces 
the path of the electron (the hole velocity being opposite 
to its momentum). The Poincare map would then con- 
sist of a single point. At non-zero excitation energy E the 
retroreflection occurs with a slight change in py , because 
of the difference 2E in the kinetic energy of electrons (at 
energy Ep + E) and holes (at energy Ep — E) . 

The resulting slow drift of the periodic trajectory 
traces out a contour in the surface of section. These 
are isochronous contours |59| . meaning that the time T 
between Andreev reflections is the same for each point 
x^Px on the contour . The adiabatic invariance of T fol- 
lows from the adiabatic invariance of the action integral 
/ over the nearly periodic motion from electron to hole 
and back to electron: 



Since _B is a constant of the motion, adiabatic invariance 
of / implies adiabatic invariance of the time T between 
Andreev reflections. 

Adiabatic invariance is defined in the limit E —^ Q 
and is therefore distinct from invariance in the sense of 
Kolmogorov-Arnold-Moser (KAM) J3, which would re- 
quire a critical E* such that a contour is exactly invari- 
ant for E < E* . There is numerical evidence 5] that 
the KAM theorem does not apply to a chaotic Andreev 
billiard. 

It is evident from Fig. ^| that contours of large T en- 
close a very small area in a chaotic system. To estimate 
the area, it is convenient to measure x in units of the 
width W of the constriction to the superconductor. Sim- 
ilarly, Px is conveniently measured in units of the range 
Ap of transverse momenta inside the constriction.^ The 
highly elongated shape evident in Fig.^lis a consequence 
of the exponential divergence in time of nearby trajecto- 
ries, characteristic of chaotic dynamics. The rate of di- 



pdq = 2ET. 



(75) 



^ We consider in this estimate the symmetric case W/L ~ 
Ap/pp <C 1, typical for the smooth confining potential of Fig. 
1151 In the asymmetric case W/L <S Ap/pp ~ 1, typical for the 
computer simulations using the kicked rot ato r, the maximal area 
Amax is smaller by a factor W/L, cf. Ref. |6r| . Consequently, the 
factor \nN in Eq. ED should be replaced by ln{NW/L). 
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FIG. 16: Poincare map for the Andreev billiard of Fig. 1151 
Each dot represents a starting point of an electron trajectory, 
at position x and with tangential momentum px (in dimen- 
sionless units). The inset shows the full surface of section, 
while the main plot is an enlargement of the central region. 
The drifting nearly periodic motion follows contours of con- 
stant time T between Andreev reflections. The cross marks 
the starting point of the trajectory shown in the previous fig- 
ure. From Ref. Issll. 



vergence is the Lyapunov exponent a. Since the Hamilto- 
nian flow is area preserving, a stretching ^_(. (i) = ^+(0)e"* 
of the dimension in one direction needs to be compen- 
sated by a squeezing £-{t) = £_(0)e~"' of the dimen- 
sion in the other direction. The area A ~ £j^£^ is then 
time independent. Initially, ^±(0) < 1. The constric- 
tion at the superconductor acts as a bottleneck, enforcing 
e±{T) < 1. These two inequalities imply £+(i) < e"(*-^), 
i- < e~"*. The enclosed area, therefore, has upper 
bound 



WApe' 



hNe-"'^, 



(76) 



where N ~ WAp/Ti 3> 1 is the number of channels in 
the point contact. 

The two invariants E and T define a two-dimensional 
torus in the four-dimensional phase space. Quantization 
of this adiabatically invariant torus proceeds following 
Einstcin-Brillouin-Keller [3, |61| , by quantizing the area 



(p pdq = 2TTh{m + ly / -i) , m^ 0,1,2, 



(77) 



enclosed by each of the two topologically independent 
contours on the torus. Eq. H77|l ensures that the wave 
functions are single valued. The integer v counts the 
number of caustics (Maslov index) and in this case should 
also include the number of Andreev reflections. 



The first contour follows the quasiperiodic orbit of Eq. 
{TSJ, leading to 



ET ^ {m + ^)TTn, m = 0,1,2,. 



(78) 



The quantization condition H78() is sufficient to determine 
the smoothed (or ensemble averaged) density of states 



{p{E)) = N / dTP{T) J2 S{E-{m^)nh/T) 



(79) 



using the classical probability distribution P{T) for the 
time between Andreev reflections. (The distribution 
P(T) is defined with a uniform measure in the surface 
of section {x,px) at the interface with the superconduc- 
tor.) 

Eq. O is the "Bohr-Sommerfeld rule" of Melsen et al. 
|2Cl|. It generalizes the familiar Bohr-Sommerfeld quan- 
tization rule for translationally invariant geometries [cf. 
Eq. Q] to arbitrary geometries. The quantization rule 
refers to classical periodic motion with period 2T and 
phase increment per period of 2ET/h— n, consisting of a 
part 2ET/Ti because of the energy difference 2E between 
electron and hole, plus a phase shift of — tt from two An- 
dreev reflections. If E is not <C A, this latter phase shift 
should be replaced by -2arccos(£'/A) [HI El IHl, cf. 
Eq. 1)21) |l. In the presence of a magnetic field an extra 
phase increment proportional to the enclosed flux should 
be included p3- Eq. H79I) can also be derived from the 
Eilenberger equation for the quasiclassical Green func- 
tion [2H . 

To find the location of individual energy levels a second 
quantization condition is needed 59]. It is provided by 
the area Srr^Pxdx enclosed by the isochronous contours, 



Pxdx ^2Trh{n + u/A), n = 0,1,2, 



(80) 



Eq. (|80|l amounts to a quantization of the period T, which 
together with Eq. H78|l leads to a quantization of E. For 
each Tn there is a ladder of Andreev levels Enm — {m + 

^)'!Th/Tn. 

While the classical T can become arbitrarily large, the 
quantized Tn has a cutoff. The cutoff follows from the 
maximal area H76|l enclosed by an isochronous contour. 
Since Eq. (|5n|l requires ^max ^ h, the longest quantized 
period is Tq = a"^[ln7V -I- 0(1)]. The lowest Andreev 
level associated with an adiabatically invariant torus is 
therefore 



Eqo = 



irh 



irha 



2Tn 21niV' 



(81) 



The time scale Tq ct I In ?i | is the Ehrenfest time te of the 
Andreev billiard, to which we will return in Sec. IVIIII 

The range of validity of adiabatic quantization is de- 
termined by the requirement that the drift 6x, Spx upon 
one iteration of the Poincare map should be small com- 
pared to the characteristic values W,pf- An estimate is 
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FIG. 17: Histograms: smoothed density of states of a bil- 
liard coupled by a ballistic A'^-mode lead to a superconductor, 
determined by Eq. II22II and averaged over a range of Fermi 
energies at fixed A'". The scattering matrix is computed nu- 
merically by matching wave functions in the billiard to trans- 
verse modes in the lead. A chaotic Sinai billiard (top inset, 
solid histogram, A^ = 20) is contrasted with an integrable 
circular billiard (bottom inset, dashed histogram, A'^ = 30). 
The solid curve is the prediction 14911 from RMT for a chaotic 
system and the dashed curve is the Bohr-Sommerfeld result 
II79I I. with dwell time distribution P{T) calculated from clas- 
sical trajectories in the circular billiard. Adapted from Ref. 
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(82) 



For low- lying levels (to ~ 1) the dimcnsionless drift is 
< 1 for r„ < Tq. Even for r„ = Tq one has 5x/W ~ 
l/lniV< 1. 



B. Integrable dynamics 

Unlike RMT, the quasiclassical theory is not restricted 
to syst ems with a chaotic classical dynamics. Melsen et 
al. J^ E3 have used the Bohr-Sommerfeld rule ((73 to 
argue that Andreev billiards with an integrable classical 
dynamics have a smoothly vanishing density of states — 
without an actual excitation gap. The presence or ab- 
sence of an excitation gap is therefore a "quantum signa- 
ture of chaos" . This is a unique property of Andreev bil- 
liards. In normal, not-superconducting billiards, it is im- 
possible to distinguish chaotic from integrable dynamics 
by looking at the density of states. One needs to measure 
density-density correlation functions for that purpose 8] . 

The difference between chaotic and integrable Andreev 
billiards is illustrated in Fig. El As expected, the chaotic 



Sinai billiard follows closely the prediction from RMT. 
(The agreement is less precise than for the kicked rotator 
of Fig. El because the number of modes iV = 20 is nec- 
essarily much smaller in this simulation.) The density of 
states of the integrable circular billiard is suppressed on 
the same mesoscopic energy scale Et as the chaotic bil- 
liard, but the suppression is smooth rather than abrupt. 
Any remaining gap is microscopic, on the scale of the 
level spacing, and therefore invisible in the smoothed 
density of states. 

That the absence of an excitation gap is generic for 
integrable billiards can be understood from the Bohr- 
Sommerfeld rule [23|. Generically, an integrable billiard 
has a power-law distribution of dwell times, P{T) ex T~p 
for r ^ oo, with p « 3 [H H^. Eq. ^ then im- 
plies a power-law density of states, {p{E)) ex E'p~^ for 
_E — > 0. The value p — i corresponds to a linearly van- 
ishing density of states. An analytical calculation ^q 
of P{T) for a rectangular billiard gives the long-time 
limit P{T) ex T^'^lnT, corresponding to the low-energy 
asymptote {p{E)) ex E\n{ET/E). The weak logarithmic 
correction to the linear density of states is consistent with 
exact quantum mechanical calculations pOl l65j| . 



C. Chaotic dynamics 

A chaotic billiard has an exponential dwell time dis- 
tribution, P{T) ex e^-^/'^'i""", instead of a power law |66| . 
(The mean dwell time is Tdweii = 2nfi/N5 = h/2ET-) 
Substitution into the Bohr-Sommerfeld rule (|79|) gives 
the density of states [69(| 



^ 2 [TrEr/Ef cosHttEt/E) 
^^^ '' 6 sinh^ (ttEt/E) 



(83) 



which vanishes ex q^'^^t/e as i? — + 0. This is a much 
more rapid decay than for integrable systems, but not 
quite the hard gap predicted by RMT ^. The two 
densities of states are compared in Fig. [TSl 

When the qualitative difference between the random- 
matrix and Bohr-Sommerfeld theories was discovered 
|20l | , it was believed to be a short-coming of the quasiclas- 
sical approximation underlying the latter theory. Lodder 
and Nazarov [21| realized that the two theoretical pre- 
dictions are actually both correct, in different limits. As 
the ratio TE/rdwcU of Ehrenfest time and dwell time is in- 
creased, the density of states crosses over from the RMT 
form H49|l to the Bohr-Sommerfeld form H83|l . We inves- 
tigate this crossover in the following section. 



VIII. QUANTUM-TO-CLASSICAL CROSSOVER 

A. Thouless versus Ehrenfest 

According to Ehrenfest's theorem, the propagation of 
a quantum mechanical wave packet is described for short 
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FIG. 18: Comparison of the smoothed density of states in 
a chaotic Andreev bilhard as it follows from RMT (Eq. (I49L 
with a hard gap) and as it follows from the Bohr-Sommerfeld 
(BS) rule (Eq. H83|l . without a hard gap). These are the two 
limiting distributions when the Ehrenfest time te is, respec- 
tively, much smaller or much larger than the mean dwell time 

Tdwell- 




FIG. 19: Two trajectories entering a chaotic billiard at a 
small separation Sx{0) diverge exponentially in time, 5x{t) — 
Sx{0)e°'^. The rate of divergence a is the Lyapunov exponent. 
An initial microscopic separation \f becomes macroscopic at 
the Ehrenfest time te = q~^ ln(L*/AF). The macroscopic 
length L* is determined by the size and shape of the billiard. 
The Ehrenfest time depends logarithmically on Planck's con- 
stant: Te = a~^ \n{Sci/h), with Sd = mvpL* the character- 
istic classical action. The evolution of a quantum mechanical 
wave packet is well described by a classical trajectory only for 
times less than te- 



times by classical equations of motion. The time scale at 
which this correspondence between quantum and classi- 
cal dynamics breaks down in a chaotic system is called 
the Ehrenfest time te JIOJ-"* As explained in Fig. ^1 
it depends logarithmically on Planck's constant: te = 
a~^ \n.{Sc\/h), with Sc\ the characteristic classical action 
of the dynamical system and a the Lyapunov exponent. 
This logarithmic /i-dependence distinguishes the 



Ehrenfest time from other characteristic time scales of 
a chaotic system, which are either /i-independent (dwell 
time, ergodic time) or algebraically dependent on h 
(Heisenberg time oc 1/(5). That the quasiclassical the- 
ory of superconductivity breaks down on time scales 
greater than te was noticed already in 1968 by Larkin 
and Ovchinnikov \I% . 

The choice of Sd depends on the physical quantity 
which one is studying. For the density of states of the An- 
dreev billiard (area A^ opening of width W <^ A^/"^ , range 
of transverse momenta Ap ~ pE inside the constriction) 
the characteristic classical action is^ Sd = mvEW^ /A^^^ 
1381. The Ehrenfest time then takes the form 



TE 



a-^[ln{N^/M) + Oil)]. 



(84) 



Here M = kEA^^^/T: and N = kFW/n are, respectively, 
the number of modes in a cross-section of the billiard 
and in the point contact. Eq. (|84|l holds for N > y/M. 
For N < vM the Ehrenfest time may be set to zero, 
because the wave packet then spreads over the entire bil- 
liard within the ergodic time |60j . 

Chaotic dynamics requires a~^ ^ TdwoU- The relative 
magnitude of te and TdwcU thus depends on whether the 
ratio N'^/M is large or small compared to the exponen- 
tially large number e"'^'*""" . 

The result of RMT 



^], cf. Sec. IvrRl is that the 
excitation gap in an Andreev billiard is of the order 
of the Thouless energy Et ~ SVrdweii- It was real- 
ized by Lodder and Nazarov |2j| that this result re- 
quires Te <C Tdwell- More generally, the excitation gap 
i^gap — min {Et, Tt-/te) is determined by the smallest of 
the Thouless and Ehrenfest energy scales. The Bohr- 
Sommerfeld theory |23, cf. Sec. lVIlO holds in the limit 
te ^ <x) and therefore produces a gapless density of 
states. 



B. Effective RMT 

A phenomenological description of the crossover from 
the Thouless to the Ehrenfest regime is p rovided by the 
"effecti ve RM T" of Silvestrov et al. ^. As described 
in Sec. I VI I Al the quasiclassical adiabatic quantization 
allows to quantize only the trajectories with periods 
T < Tq = Te- The excitation gap of the Andreev billiard 
is determined by the part of phase space with periods 
longer than te- Effective RMT is based on the hypoth- 
esis that this part of phase space can be quantized by a 
scattering matrix ScS in the circular ensemble of RMT, 
with a reduced dimensionality 



N, 



cff 



N 



I P{T) dT = ATe-^^/^d, 



(85) 



** The name "Ehrenfest time" was coined in Ref. I7ll . 



^ The simpler expression Sd = thveW of Ref. l5fll applies to the 
symmetric case W/A^''^ ~ Ap/pp <SC 1. 
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FIG. 20: Pictorial representation of the effective RMT of an 
Andreev billiard. The part of phase space with time T > te 
between Andreev reflections is represented by a chaotic cavity 
(mean level spacing (5cff), connected to the superconductor by 
a long lead {N^s propagating modes, one-way delay time te/2 
for each mode). Between two Andreev reflections an electron 
or hole spends, on average, a time te in the lead and a time 
Tdwcii in the cavity. The scattering matrix of lead plus cav- 
ity is exp{iETE /h)So(E), with So{E) distributed according 
to the circular ensemble of RMT (with effective parameters 
A'^eff , 5oft). The complete excitation spectrum of the Andreev 
billiard consists of the levels of the effective RMT (periods 
> Te) plus the levels obtained by adiabatic quantization (pe- 
riods < Te). 



The energy dependence of Scs{E) is that of a chaotic 
cavity with mean level spacing Scs, coupled to the super- 
conductor by a long lead with Nog propagating modes. 
(See Fig. 1201) The lead introduces a mode-independent 
delay time te between Andreev reflections, to ensure that 
P{T) is cut off for T < te- Because P{T) is exponential 
ex exp(— T/rdwoii), the mean time (T)* between Andreev 
reflections in the accessible part of phase space is sim- 
ply Te + Tdwcii ■ The effective level spacing in the chaotic 
cavity by itself (without the lead) is then determined by 



2Trn 



off 



= (T)* -Te = Td- 



well- 



(86) 



It is convenient to separate the energy dependence 
due to the lead from that due to the cavity, by writ- 
ing SesiE) = eyip{iETE/^)So{E), where So{E) repre- 
sents only the cavity and has an energy dependence of 
the usual RMT form H30(l — with effective parameters 
iVofj and Jeff- The determinant equation (|22|l for the ex- 
citation spectrum then takes the form 

Dot \l - a{Efe^''^^^/'''So{E)Soi-Ey] =0. (87) 

We can safely replace a{E) = exp[— i arccos(i?/A)] -^ —i 
(since E <^ A), but the energy dependence of the phase 
factor e'^^^'^E/fi ^g^y^ j^q^ ^g omitted. 

In App. we calculate the smallest positive E that 
solves Eq. (|87|l , which is the excitation gap -Egap of the ef- 
fective RMT. The result is plotted in Fig.|^(sohd curve), 
as a function of r^j/rdwoii- The two asymptotes (dotted 
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FIG. 21: Excitation gap of the Andreev billiard in the 
crossover from Thouless to Ehrenfest regimes. The solid curve 
is the solution of the effective RMT of Ref. ^j, derived in 
App. ^ The dotted lines are the two asymptotes 1881 and 
18911 . The dashed curve is the result of the stochastic model 
of Ref. H, discussed in Sec lVniOl 
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2^ 



, Te ^ Tdwcl(?8) 



l-(3 + V8)^^^V TE» 
Te ) 



T'dwcu,(89) 



with 7 = ^(vS — 1) the golden number. 

The Te time delay characteristic of the effective RMT 
was introduced in Ref. [53 , but its effect on the excitation 
gap was not evaluated properly.^ As a consequence the 
formula for the gap given in that paper, 



E, 






0.30h 



TE + Tdwell 



(90) 



provides only a qualitative description of the actual 
crossover. 

The inverse correlation (|9()|l between the gap and the 
dwell time of long trajectories was observed in a computer 
simulation of the Andreev kicked rotator ^3 ■ The data 
points in Fig. 1221 track the excitation gap as the location 
in phase space of the NS interface is varied. The solid 
curve is a plot of 



1 _ j^.P{T)dT 



(T)* J^,TP{T)dT' 



(91) 



I am indebted to P. W. Brouwer for spotting the error. 
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FIG. 22: The data points (left axis) are the quantum me- 
chanical gap values -Egap of the Andreev kicked rotator as a 
function of the location of the NS interface, for parameter val- 
ues M = 131072, Tdwoii = M/N = 5, K = 14. The solid curve 
(right axis) is the reciprocal of the mean dwell time 19111 of 
classical trajectories longer than T* — 7. Adapted from Ref. 



with P{T) the classical dvifcU time distribution and T* = 
7. We see that the sample-to-sample fluctuations in the 
gap correlate very well with the fluctuations in the mean 
dwell time of long trajectories. The correlation is not 
sensitive to the choice of T* , as long as it is greater than 
TE = 4.4. 



C. Stochastic model 

Small-angle scattering by a smooth disorder potential 
provides a stochastic model for the quantum diffraction 
of a wave packet in a chaotic billiard 74] . The scattering 
time of the stochastic model plays the role of the Ehren- 
fest time in the deterministic chaotic dynamics. The ad- 
vantage of a stochastic description is that one can average 
over different realizations of the disorder potential. This 
provides for an established set of analytical techniques. 
The disadvantage is that one does not know how well 
stochastic scattering mimics quantum diffraction. 

Vavilov and Larkin 38] have used the stochastic model 
to study the crossover from the Thouless regime to the 
Ehrenfest regime in an Andreev billiard. They discov- 
ered that the rapid turn-on of quantum diffraction at 
te <] Tdwcii not only causes an excitation gap to open 
at Ti/te, but that it also causes oscillations with period 
Ti/te in the ensemble-averaged density of states {p{E)) 
at high energies E > Ex- In normal billiards oscillations 
with this periodicity appear in the level-level correlation 




FIG. 23: Oscillatory density of states at finite Ehrenfest 
time (solid curve), compared with the smooth limits of zero 
(RMT) and infinite (BS) Ehrenfest times. The solid curve is 
the result of the stochastic model of Vavilov and Larkin, for 
Te — 'i Tdwcii = 2iTi/2Et. (The definition l|84|l of the Ehrenfest 
time used here differs by a factor of two from that used by 
those authors.) Adapted from Ref. |3q[. 



function |73, but not in the level density itself. 

The predicted oscillatory high-energy tail of {p{E)) is 
plotted in Fig. [531 for the case r£;/rdweii = 3, together 
with the smooth results of RMT {te /rdwcW -^ 0) and 
Bohr-Sommerfeld (BS) theory (r£:/Tdweii —* oo). 

Independent analytical support for the existence of os- 
cillations in the density of states with period TT'/te comes 
from the singular perturbation theory of Ref. 76] . Sup- 
port from numerical simulations is still lacking. Jacquod 
et al. [23 did find pronounced oscillations for E > Et in 
the level density of the Andreev kicked rotator. However, 
since these could be described by the Bohr-Sommerfeld 
theory they can not be the result of quantum diffraction, 
but must be due to nonergodic trajectories |77j . 

The T£;-dependence of the gap obtained by Vavilov and 
Larkin is plotted in Fig. [2] (dashed curve). It is close 
to the result of the effective RMT (solid curve). The 
two theories predict the same limit -Egap —>■ ttTjI'Itk for 
TB/Tdwcii -^ 00. The asymptotes given in ReL |33] are 



E. 



E, 



TdwcU 



, TE < rd„oii,(92) 



2te \ TE 



2TdwcU , 

TE > Tdwoll- (93) 



Both are different from the results (|88|) and (|89|l of the 
effective RMT.'^ 



^ Since 27 — 1 = 0.236, the small- r^ asymptote of Vavilov and 
Larkin differs by a factor of two from that of the effective RMT. 
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FIG. 24: Ehrenfest-time dependence of the excitation gap in 
an Andreev billiard, according to the effective RMT (solid 
curve, calculated in App.^ and according to the stochastic 
model (dashed curve, calculated in Ref. 38]). The data points 
result from the simulation of the Andreev kicked rotator |25|| 
(closed circles, in the range A'^ = 10^ — 10^) and of the Sinai 
billiard shown in the inset |7q| (open circles, in the range 
A' = 18 - 30). 



D. Numerical simulations 

Because the Ehrenfest time grows only logarithmically 
with the size of the system, it is exceedingly difficult to 
do numerical simulations deep in the Ehrenfest regime. 
Two simulations |2^ [Tg have been able to probe the 
initial decay of the excitation gap, when te ^ TdwoU- We 
show the results of both simulations in Fig. |31 (closed and 
open circles), together with the full decay as predicted by 
the effective RMT of Se c.lVllllJI fsolid curve) and by the 
stochastic model of Sec. IVIII CI (dashed curve). 

The closed circles were obtained by Jacquod et al. |23 
using the stroboscopic model of Sec. (the Andreev 
kicked rotator). The number of modes N in the con- 
tact to the superconductor was increased from 10^ to 
10^ at fixed dwell time TdwcU = M/N = 5 and kicking 
strength K = U (corresponding to a Lyapunov exponent 
a « ln(i^/2) = 1.95). In this way all classical properties 
of the billiard remain the same while the effective Planck 
constant hcs — 1/M = 1/A^rd„cii is reduced by three 
orders of magnitude. To plot the data as a function of 
TB/Tdwoii, Eq- (|84|l was used for the Ehrenfest time. The 
unspecified terms of order unity in that equation were 
treated as a single fit parameter. (This amounts to a 
horizontal shift by —0.286 of the data points in Fig. 1241) 

The open circles were obtained by Kormanyos et al. 
|j8] for the chaotic Sinai billiard shown in the inset. The 
number of modes N was varied from 18 to 30 by varying 
the width of the contact to the superconductor. The 
Lyapunov exponent a ~ 1.7 was fixed, but TdwoU was 
not kept constant in this simulation. The Ehrenfest time 
was computed by means of the same formula (|84|l . with 
M = 2Lckp/TT and Lc the average length of a trajectory 



between two consecutive bounces at the curved boundary 
segment. 

The data points from both simulations have substan- 
tial error bars (up to 10%). Because of that and because 
of their limited range, we can not conclude that the sim- 
ulations clearly favor one theory over the other. 



IX. CONCLUSION 

Looking back at what we have learned from the study 
of Andreev billiards, we would single out the breakdown 
of random-matrix theory as the most unexpected discov- 
ery and the one with the most far-reaching implications 
for the field of quantum chaos. In an isolated chaotic bil- 
liard RMT provides an accurate description of the spec- 
tral statistics on energy scales below h/Tc^g (the inverse 
ergodic time). The weak coupling to a superconductor 
causes RMT to fail at a much smaller energy scale of 
^/Td-weii (the inverse of the mean time between Andreev 
reflections), once the Ehrenfest time te becomes greater 
than TdwcU. 

In the limit te — > oo, the quasiclassical Bohr- 
Sommerfeld theory takes over from RMT. While in iso- 
lated billiards such an approach can only be used for 
integrable dynamics, the Bohr-Sommerfeld theory of An- 
dreev billiards applies regardless of whether the classical 
motion is integrable or chaotic. This is a demonstration 
of how the time-reversing property of Andreev reflection 
unravels chaotic dynamics. 

What is lacking is a conclusive theory for finite te ^ 
TdwcU- The two phenomenological approaches of Sees. 
IVIII Bl and IVIII HI agree on the asymptotic behavior 



lim E'^an = 



nha 



h-*o 



gap 



2| ln?i| -I- constant' 



(94) 



in the classical ?i ^ limit (understood as iV ^ oo at 
fixed TdwcU ). There is still some disagreement on how 
this limit is approached. We would hope that a fully 
microsco pic approach, for example based on the ballistic 
(T- model |Z^|s3) could provide a conclusive answer. At 
present technical difficulties still stand in the way of a 
solution along those lines (Slj . 

A new direction of research is to investigate the effects 
of a nonisotropic superconducting order parameter on the 
Andreev billiard. The case of d-wave symmetry is most 
interesting because of its relevance for high-temperature 
superconductors. The key ingredients needed for a theo- 
retical description exist, notably RMT |S3|, quasiclassics 
|83l |. and a numerically efficient Andreev map |84l |. 
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still leads to Eq. (|43a|) . but Eq. 
by 



I43b|) should be replaced 



APPENDIX A: EXCITATION GAP IN 

EFFECTIVE RMT AND RELATIONSHIP WITH 

DELAY TIMES 

We seek the edge of the excitation spectrum as it fol- 
lows from the determinant equation H87|) . which in zero 
magnetic field and for E <^ A takes the form 



Det 



1 



2iE 



^^^^So{E)So{~Ey 



0. 



(Al) 



The unitary symmetric matrix ^o has the RMT distri- 
bution of a chaotic cavity with effective parameters A^cff 
and i5eff given by Eqs. ifH^ and (IHEJ- The mean dwell 
time associated with Sq is TdwoU- The calculation for 
Ncs ^ 1 follows the method described in Sees. IVI^ and 
modified as in Ref. 



IVIBI modified as in Ref. |35j to account for the energy 
dependent phase factor in the determinant. 

Since 5*0 is of the RMT form H3U|) . we can write Eq. 
(jAip in the Hamiltonian form (|32|l . The extra phase 
factor exjp(2iETE/fi) introduces an energy dependence of 
the coupling matrix, 



WoiE) 



cosu 



WqW^ sin u 



WoW^sini 



(A2) 



where we have abbreviated u = Ete/TT'- The subscript 
reminds us that the coupling matrix refers to the reduced 
set of iVeff channels in the effective RMT. Since there is no 
tunnel barrier in this case, the matrix Wq is determined 
by Eq. (EU with r„ = 1. The Hamiltonian 



TLo 




(A3) 



is that of a chaotic cavity with mean level spacing S^fi- 
We seek the gap in the density of states 

p{E) = -^ImTr M + ^ j {E + zO+ - T^o + >Vo)-\ 

(A4) 
cf. Eq. (p|l. 

The selfconsistency equation for the ensemble-averaged 
Green function, 

g=[E + Wo- {MS,s/7r)a,Ga,]-\ (A5) 



Gll+Gi2SinM = -(Tdwcll/TB)MGi2 

X (Gi2 + cosu + Gil sinu).(A6) 

(We have used that NosScS = 27r?i/Tdweii-) Because of the 
energy dependence of the coupling matrix, the equation 
(|H|l for the ensemble averaged density of states should 
be replaced by 



(p(i?))---Im(Gii 



COSM 



-Gi 



(A7) 



The excitation gap corresponds to a square root sin- 
gularity in (p(_E)), which can be obtained by solving 
Eqs. H43a() and IjAGp jointly with dE/dGn = for 
u £ (0,7r/2). The result is plotted in Fig. EH The small- 
and large- tb asymptotes are given by Eqs. (|SS|l and l|S!?|l . 

The large- r^; asymptote is determined by the largest 
eigenvalue of the time-delay matrix. To see this rela- 
tionship, note that for te ^ Tdwoii we may replace the 
determinant equation ljAl|l by 



Dot 



1 + cxp[2iETE/n + 2iEQ{Q)] + O 



T'dwell 
TE 



The Hermitian matrix 



Q{E)^-So{E)^^So{E) 
I aE 



= 0. 

(A8) 

(A9) 



is known in RMT as the Wigner-Smith or time-delay ma- 
trix. The roots Enp of Eq. IjASp satisfy 



2Enp{TE + r„) = (2p + l)7rn, p = 0, 1, 2, , 



(AlO) 



The eigenvalues ti , r2 , . . . T^^f^ of hQ{0) are the delay 
times. They are all positive, distributed according to a 
generalized Laguerre ensemble [S^. In the limit A'cir -^ 
oo the distribution of the Tn 's is nonzero only in the inter- 
val (t_, T_|_), with T± = T(iweii(3 ± v^). By taking p = 0, 
Tn — T+ we arrive at the asymptote (|89|l . 

The precise one-to-one correspondence (JAIOII between 
the spectrum of low-lying energy levels of the Andreev 
billiard and the spectrum of delay times is a special prop- 
erty of the classical limit r^ ^ oo. For te ^ Tdwcii there 
is only a qualitative similarity of the two spectral densi- 
ties ISGl. 



[1] A. F. Andreev, Zh. Eksp. Teor. Fiz. 46, 1823 (1964) [Sov. 

Phys. JETP 19, 1228 (1964)]. 
[2] Y. Imry, Introduction to Mesoscopic Physics (Oxford 

University, Oxford, 2002). 
[3] C. W. J. Beenakker, Rev. Mod. Phys. 69, 731 (1997). 



[4] B. J. van Wees and H. Takayanagi, in Mesoscopic Elec- 
tron Transport, edited by L. L. Sohn, L. P. Kouwenhoven, 
and G. Schon, NATO ASI Series E345 (Kluwer, Dor- 
drecht, 1997). 

[5] I. Kosztin, D. L. Maslov, and P. M. Goldbart, Phys. Rev. 



21 



Lett. 75, 1735 (1995). 
[6] J. Eroms, M. Tolkiehn, D. Weiss, U. Rossler, J. De Boeck, 

and G. Borghs, Europhys. Lett. 58, 569 (2002). 
[7] M. C. Gutzwiller, Chaos in Classical and Quantum Me- 
chanics (Springer, New York, 1990). 
[8] F. Haake, Quantum Signatures of Chaos (Springer, 

Berlin, 2001). 
[9] M. L. Mehta, Random Matrices (Academic, New York, 

1991). 
[10] T. Guhr, A. Miiller-Groeling, and H. A. Weidenmiiller, 

Phys. Rep. 299, 189 (1998). 
[11] P. G. de Gennes, Superconductivity of Metals and Alloys 

(Benjamin, New York, 1966). 
[12] K. K. Likharev, Rev. Mod. Phys. 51, 101 (1979). 
[13] C. W. J. Beenakker, in: Transport Phenomena in Meso- 

scopic Systems, edited by H. Fukuyama and T. Ando 

(Springer, Berlin, 1992); cond-mat/0406127 
[14] W. L. McMillan, Phys. Rev. 175, 537 (1968). 
[15] P. G. de Gennes and D. Saint-James, Phys. Lett. 4, 151 

(1963). 
[16] A. A. Golubov and M. Yu. Kuprianov, J. Low Temp. 

Phys. 70, 83 (1988). 
[17] W. Belzig, C. Bruder, and G. Schon, Phys. Rev. B 54, 

9443 (1996). 
[18] S. Pilgram, W. Belzig, and C. Bruder, Phys. Rev. B 62, 

12462 (2000). 
[19] W. Belzig, F. K. Wilhelm, C. Bruder, G. Schon, and A. 

D. Zaikin, Superlatt. Microstruct. 25, 1251 (1999). 
[20] J. A. Melsen, P. W. Brouwer, K. M. Frahm, and C. W. 

J. Beenakker, Europhys. Lett. 35, 7 (1996). 
[21] A. Lodder and Yu. V. Nazarov, Phys. Rev. B 58, 5783 

(1998). 
[22] C. W. J. Beenakker, Phys. Rev. Lett. 67, 3836 (1991); 

68, 1442(E) (1992). 
[23] E. B. Bogomolny, Nonlinearity 5, 805 (1992). 
[24] R. E. Prange, Phys. Rev. Lett. 90, 070401 (2003). 
[25] Ph. Jacquod, H. Schomerus, and C. W. J. Beenakker, 

Phys. Rev. Lett. 90, 207004 (2003). 
[26] A. Ossipov, T. Kottos, and T. Geisel, Europhys. Lett. 

62, 719 (2003). 
[27] Y. V. Fyodorov and H.-J. Sommers, JETP Lett. 72, 422 

(2000). 
[28] A. Ossipov and T. Kottos, Phys. Rev. Lett. 92, 017004 

(2004). 
[29] F. M. Izrailev, Phys. Rep. 196, 299 (1990). 
[30] J. Tworzydlo, A. Tajic, and C. W. J. Beenakker, 

'cond-mat/0405122' 
[31] R. Ketzmerick, K. Kruse, and T. Geisel, Physica D 131, 

247 (1999). 
[32] O. Bohigas, M.-J. Giannoni, and C. Schmit, Phys. Rev. 

Lett. 52, 1 (1984). 
[33] S. Miiller, S. Heusler, P. Braun, F. Haake, and A. Ah- 

land, nlin. CD/0401021 
[34] M. G. Vavilov, P. W. Brouwer, V. Ambegaokar, and C. 

W. J. Beenakker, Phys. Rev. Lett. 86, 874 (2001). 
[35] P. W. Brouwer and C. W. J. Beenakker, Chaos, Solitons 

& Fractals 8, 1249 (1997). In Eq. (8) of this publication 

the factor 2kT should read 4kT and in the right-hand- 
side of Eq. (15) a factor 1 -I- dW /de should be inserted 

(cf. Eq. 1351 in the present paper). 
[36] K. M. Frahm, P. W. Brouwer, J. A. Melsen, and C. W. 

J. Beenakker, Phys. Rev. Lett. 76, 2981 (1996). 
[37] A. Ahland and M. R. Zirnbauer, Phys. Rev. Lett. 76, 

3420 (1996); Phys. Rev. B 55, 1142 (1997). 



[38] M. G. Vavilov and A. I. Larkin, Phys. Rev. B 67, 115335 

(2003). 
[39] A. Pandey, Ann. Phys. (N.Y.) 134, 110 (1981). 
[40] E. Brezin and A. Zee, Phys. Rev. E 49, 2588 (1994). 
[41] L. A. Pastur, Teoret. Mat. Fiz. 10, 102 (1972) [Theoret. 

Math. Phys. 10, 67 (1972)]. 
[42] M. Tinkham, Introduction to Superconductivity 

(McGraw-HiU, New York, 1995). 
[43] J. A. Melsen, P. W. Brouwer, K. M. Frahm, and C. W. 

J. Beenakker, Physica Scripta T69, 223 (1997). 
[44] A. Pandey and M. L. Mehta, Commun. Math. Phys. 87, 

449 (1983). 
[45] T. Nagao and K. Slevin, J. Math. Phys. 34, 2075 (1993). 
[46] J. T. Bruun, S. N. Evangelou, and C. J. Lambert, J. 

Phys. Condens. Matt. 7, 4033 (1995). 
[47] K. Efetov, Supersymmetry in Disorder and Chaos (Cam- 
bridge University, Cambridge, 1997). 
[48] S. Gnutzmann, B. Self, F. von Oppen, and M. R. Zirn- 
bauer, Phys. Rev. E 67, 046225 (2003). 
[49] C. W. J. Beenakker, Phys. Rev. B 46, 12841 (1992). 
[50] C. A. Tracy and H. Widom, Commun. Math. Phys. 159, 

151 (1994); 177, 727 (1996). 
[51] P. M. Ostrovsky, M. A. Skvortsov, and M. V. Feigelman, 

Phys. Rev. Lett. 87, 027002 (2001); JETP Lett. 75, 336 

(2002). 
[52] A. Lamacraft and B. D. Simons, Phys. Rev. B 64, 014514 

(2001). 
[53] P. M. Ostrovksy, M. A. Skvortsov, and M. V. Feigelman, 

Phys. Rev. Lett. 92, 176805 (2004). 
[54] L. G. Aslamasov, A. L Larkin, and Yu. N. Ovchinnikov, 

Zh. Eksp. Teor. Fiz. 55, 323 (1968) [Sov. Phys. JETP 

28, 171 (1969)]. 
[55] A. V. Shytov, P. A. Lee, and L. S. Levitov, Phys. Uspekhi 

41, 207 (1998). 
[56] J. Wiersig, Phys. Rev. E 65, 036221 (2002). 
[57] I. Adagideh and P. M. Goldbart, Phys. Rev. B 65, 201306 

(2002). 
[58] J. Cserti, P. Polinak, G. Palla, U. Ziilicke, and C. J. 

Lambert, Phys. Rev. B 69, 134514 (2004). 
[59] P. G. Silvestrov, M. C. Goorden, and C. W. J. Beenakker, 

Phys. Rev. Lett. 90, 116801 (2003). 
[60] P. G. Silvestrov, M. C. Goorden, and G. W. J. Beenakker, 

Phys. Rev. B 67, 241301 (2003). 
[61] K. P. Duncan and B. L. Gyorffy, Ann. Phys. (N.Y.) 298, 

273 (2002). 
[62] J. Cserti, A. Kormanyos, Z. Kaufmann, J Koltai, and C. 

J. Lambert, Phys. Rev. Lett. 89, 057001 (2002). 
[63] J. Cserti, A. Bodor, J. Koltai, and G. Vattay, Phys. Rev. 

B 66, 064528 (2002). 
[64] J. Cserti, B. Beri, P. PoUner, and Z. Kaufmann, 

cond-mat/0405404 
[65] W. Ihra, M. Leadbeater, J. L. Vega, and K. Richter, Eu- 
rophys. J. B 21, 425 (2001). 
[66] W. Bauer and G. F. Bertsch, Phys. Rev. Lett. 65, 2213 

(1990). 
[67] H. U. Baranger, R. A. Jalabert, and A. D. Stone, Chaos 

3, 665 (1993). 
[68] A. Kormanyos, Z. Kaufmann, J. Cserti, and C. J. Lam- 
bert, Phys. Rev. B 67, 172506 (2003). 
[69] H. Schomerus and C. W. J. Beenakker, Phys. Rev. Lett. 

82, 2951 (1999). 
[70] G. P. Berman and G. M. Zaslavsky, Physica A 91, 450 

(1978). 
[71] B. V. Chirikov, F. M. Izrailev, and D. L. Shepelyansky, 



22 



Physica D 33, 77 (1988). [80] 

[72] A. I. Larkin and Yu. N. Ovchinnikov, Zh. Eksp. Teor. 

Fiz. 55, 2262 (1968) [Sov. Phys. JETP 28, 1200 (1969)]. [81] 

[73] M. C. Goorden, Ph. Jacquod, and C. W. J. Beenakker, 

Phys. Rev. B 68, 220501 (2003). [82] 

[74] I. L. Aleiner and A. I. Larkin, Phys. Rev. B 54, 14423 

(1996). [83] 

[75] I. L. Aleiner and A. I. Larkin, Phys. Rev. E 55, 1243 

(1997). [84] 

[76] L Adagideh and C. W. J. Beenakker, Phys. Rev. Lett. 

89, 237002 (2002). [85] 

[77] W. Ihra and K. Richter, Physica E 9, 362 (2001). 
[78] A. Kormanyos, Z. Kaufmann, C. J. Lambert, and J. 

Cserti, cond-mat/0309306 [86] 

[79] B. A. Muzykantskii and D. E. Khmelnitskii, JETP Lett. 

62, 76 (1995). 



A. V. Andreev, B. D. Simons, O. Agam, and B. L. Alt- 
shuler, Nucl. Phys. B 482, 536 (1996). 
D. Taras-Semchuk and A. Altland, Phys. Rev. B 64, 
014512 (2001). 

A. Altland, B. D. Simons, and M. R. Zirnbauer, Phys. 
Rep. 359, 283 (2002). 

I. Adagideh, D. E. Sheehy, and P. M. Goldbart, Phys. 
Rev. B 66, 140512 (2002). 

I. Adagideh and Ph. Jacquod, Phys. Rev. B 69, 020503 
(2004). 

P. W. Brouwer, K. M. Frahm, and C. W. J. Beenakker, 
Phys. Rev. Lett. 78, 4737 (1997); Waves in Random Me- 
dia, 9, 91 (1999). 

M. G. A. Crawford and P. W. Brouwer, Phys. Rev. E 65, 
026221 (2002). 



